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1 Introduction 



At the dawn of quantum mechanics physicists realized that not every partial 
differential equation nature has to offer can be solved analytically. Even the 
most fundamental equation of quantum mechanics, the Schrodinger equation, 
can be solved analytically only for very few systems. As Dirac put it in 1929, 
"the underlying physical laws necessary for the mathematical theory of a 
large part of physics and the whole of chemistry are thus completely known, 
and the difficulty is only that the exact application of these laws leads to 
equations much too complicated to be soluble". Thus, to apply quantum 
mechanical concepts to generic molecular systems one has to find some other 
way of solving the complicated equations. 

One approach is to solve differential equation numerically, via introduction 
of a finite grid. This reduces the problem to having to solve a (large) system 
of linear equations, a problem well suited for solving on digital computers. 
However, high accuracy via this approach can be achieved only with non- 
uniform grids, with high density of points around the nuclei. 

An alternative approach, more intuitively obvious for chemists familiar with 
the concept of atomic orbitals (physicists familiar with Bloch functions), is 
to expand the solution in terms of functions that behave like the solution, 
atomic orbitals in the chemist's case (plane waves in the physicist's). Again, 
this simplifies the problem down to solving matrix equations, with one com- 
plication being that the matrix elements become more or less non-trivial to 
compute. The matrix elements of the Hamiltonian operator consist of several 
types of molecular integrals, and the evaluation of these integrals is the main 
focus of this document. 

Of course, one has to anticipate what the exact solution looks like, but it turns 
out that the fundamental properties of atomic and molecular wavefunctions, 
such as nuclear and electron-electron cusps, exponential decay at infinity, 
etc. are known and it is relatively easy to find trial functions to satisfy 
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them. Atom-centered basis sets help to describe the critical region around 
the nuclei efficiently, hence these have been used in a majority of quantum 
chemical studies of molecules. (Homework problem - find a study that didn't 
utilize atom-centered functions...) 

The choice of a particular functional form for the expansion is more difficult. 
If one had to solve the Schrodinger equation for the hydrogen atom in finite 
basis, the natural choice would be to choose functions of the form 



that resemble the well-known discrete-spectrum solutions found in any text- 
book on quantum mechanics. Introduced by Slater, this type of basis func- 
tions is particularly well suited for atomic calculations. However, the ap- 
plicability Slater-type functions to other than trivial molecular problems is 
hindered by the enormous computational complexity of the resulting expres- 
sions for matrix elements of the Hamiltonian. 

Boys found a more tractable choice in 1950 1 as he introduced Gaussian-type 
functions of the form 



and they have been with us ever since. Much simpler expressions for the 
matrix elements more than compensate for improper behavior of Gaussians 
at the origin and infinity. Indeed, an s-type Gaussian (a Gaussian functions 
with / + m + n equal 0) is smooth at the origin, whereas an s-type Slater- type 
function has a cusp at the origin (non-zero derivative with respect to r). Also, 
Gaussian-type functions decay with r much faster that Slater-type functions. 
However, a fact that a given Slater-type functions can be well represented 
as a linear combination of only few Gaussians with different exponents was 
noticed early on, thus STO-nG basis sets were introduced, in which a single 
Slater-type function is represented as a linear combination of n Gaussians. 

One important note to make before we move on to mathematical details is 

1 S. F. Boys, Proc R. Soc. London Ser. A 200, 542 (1950). 



= r L e-< r Yl'(6, <f) 



(1.1) 



(j) = x l y m z n e 
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the value of efficient algorithms for computing integrals. Normally, in a typ- 
ical high accuracy calculation only a small portion of CPU time is spent in 
computing molecular integrals, and the major part is spent in computing 
wavefunction parameters. However, the situation has changed dramatically 
since the late 1980s. Direct methods which do not store the integrals but re- 
compute them as needed put the integrals evaluation process into a spotlight 
again. OK, enough said, on to the integrals! 

2 Elementary Basis Function Analysis 
2.1 Normalization 

As we have mentioned, the standard basis functions used in ab initio theory 
are Gaussian functions, or linear combinations thereof. There are two types 
of Gaussians: Cartesian and spherical harmonic. The functional expression 
for the unnormalized primitive Cartesian Gaussian-type functions is worth 
rewriting here: 

M (r) = x l y m z n e- ar2 (2.1) 

This is a basis function of angular momentum (I + m + n), 2 centered at the 
origin, with orbital exponent a. The term "primitive" denotes a fact that it 
is a single function. A linear combination of primitive Gaussians located at 
the same center 

c 

M ( r) = ^ x h y mi z n 'e~ a ^ (2.2) 

i=0 

is called a contracted Cartesian Gaussian function. Standard notation. A cor- 
responding primitive spherical harmonic, or pure angular momentum, Gaus- 
sian is simply 

^(r) = rV-V L M (^) (2.3) 

Note that we use capital letters for the angular momentum number and 
the projection of angular momentum on z-axis to distinguish those from 

2 Angular momentum is a slippery term here, some authors just use "orbital quantum number" 
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the exponents of x and y in Eqn (2.1). Hence, L is related to / + m + 
n, but not identical. The rule is that any spherical harmonic Gaussian of 
angular momentum L can be expressed solely in terms Cartesian Gaussians of 
angular momentum L. However, the reverse is not generally true. Similarities 
end there though. Just keep in mind that unless stated otherwise, terms 
"Gaussian" and "Gaussian function" will refer to Cartesian Gaussian-type 
functions throughout this document. 

These atomic orbital-like basis functions need not be orthogonal to one an- 
other, but for later convenience, it would be nice to have them normalized. 
Thus impose the condition 

J 0;(r)^(r)dr = 1. (2.4) 

Let's evaluate this integral. Assume a normalization constant of N for M , 
and call (2.4) a self-overlap integral, SO. 

SO = J N 2 {x l y m z n e- ar2 ){x l y m z n e- ar2 )dr (2.5) 

This integral over all space is separable when done in Cartesian coordinates. 3 
Using r 2 = x 2 + y 2 + z 2 and dr = dxdydz, we get 

SO = N 2 J dxx 2l e- 2ax2 J dyy 2m e- 2ay2 J dzz 2n e- 2az2 (2.7) 

= N 2 IJ y I z (2.8) 

Full derivation of these integrals (I x , etc) is left to the reader (otherwise this 
would not be a learning experience). The result is 



oo 



21 -2a X > (2/-l)!!v^ 



I x = / dxx M e- lax = ± (2.9) 



Spherical harmonic Gaussians aren't as nice in this respect, that's why normally all integrals are computed in 
the Cartesian basis and then transformed into the spherical harmonic basis. 
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= 1 



Recall that {21 - 1)!! = 1 • 3 • 5 • • • (2/ - 1). Thus 

= 2 r (2/-l)!!(2m-l)!!(2n-l)!!7r 3 / 2 

Rearranging that to solve for N, the normalization constant, 

^2\ 3 / 4 2( /+m + n )o/( 2/+2m+2n+3 )/ 4 
V = [n) [(2/-l)!!(2m-l)!!(2n-l)!!] 1 /2 



(2.10) 



(2.11) 



This result is completely general - for uncontracted functions. As you might 
have guessed, computing normalization constants for contracted Gaussians 
is not much more difficult. 



2.2 Products of Contracted Cartesian Gaussians 



Examine some contracted s functions. Let 



n 



0(r) = Nj2 a * e ~ a * r " 



(2.12) 



where n is the number of primitive functions in the contracted function 0(r), 
and cii are the contraction coefficients. The product 0*(r)0(r) can be written 



0*(r)0(r) = N z 



n 



a 3 e 



—a^r 



(2.13) 



Since the bracketed term contains a product of two polynomials, only two 
types of terms can result; the square of each uncontracted function and the 
biproducts of different uncontracted functions. Take an example where n is 
three: 

[aie- a ^ 2 + a 2 e- a2r2 + a 3 e^ r2 } 2 = [a\e- 2a ^ + a 2 ^ 2 ^ 2 + a\e- 2a ^ 

+2a 2 a 3 e- {a " +a3)r2 } (2.14) 
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In this case, as in all others, there are only two types of terms of which the 
integral needs to be taken. They may be written and evaluated as 



1. aje- 2a * r2 dr = aj 



IT 



2a, 



3/2 



2. / 2a l a J e~ {ai+ ^ )r2 dv = 2 ai a 3 



IT 



3/2 



a,- + a 



(2.15) 
(2.16) 



It is realized that 1. above can be obtained by setting i = j in 2., and 
henceforth only the general case needs to be considered. Generalizing to ar- 
bitrary n is straightforward, and so the normalization of contracted Gaussian 
functions can proceed as 



J 0*(r)0(r)dr = iVV /2 



(2c*i) 3 / 2 

n n 



+ ••• + 



CliClj 



[a x + a 2 ) 3 / 2 
= 1, 



+ 



^ (a t + a 7 -) 3/2 
j 

thus the normalization constant for the entire contraction will be 

-1/2 



(2.17) 



N = 7T- 3 / 4 



E 



(2.18) 



Contractions of Gaussians of arbitrary angular momentum are a bit worse, 
but if we assume all of the contracted functions to be of the same angular 
momentum 
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0(r) = N 



a lX l y m z n e^ r2 + a 2 x l y m z n e- a ^ + 



= Nx l y m z n J2 a i e 



J 0*(r)0(r)dr = N 2 J x 2l y 2m z 
J 0*(r)0(r)dr 



^ 2 EE 

i=0 j=0 



diCLj i \x~y 



x 2l 2m z 2n e -a^ e -a^ 



(2.19) 
dr(2.20) 
dr (2.21) 



The product in brackets in Eqn (2.21) we've encountered before. Analogous 
to Eqn (2.9), the general form for the integral in the double sum is 



x U y 2rn z 2n e -(a i+aj y dr = ^ 3/2 



(2/-l)!!(2m-l)!!(2n-l)!! 



(2.22) 



The self overlap is then 

N 2 n 3 / 2 (2l - l)!!(2m - l)!!(2n - 1)!! 



0*(r)0(r)dr = 



2l+m+n 



E 



i,3 



CliClj 



(a l + a j y+ m+n +' i / 2 ' 



(2.23) 

Calling / + m + n = L, the angular momentum of the shell, and solving for 



^ _ N 2 ^' 2 {21 - l)!!(2m - l)!!(2n - 1)!! ^ 



didj 



2 L 



h3 



i (a t + a,) L+3 / 2 



= 1 



N = 



7r 3 / 2 (2/-l)!!(2m-l)!!(2n-l)!! 



E 



1-1/2 



(2.24) 
(2.25) 



Before we approach one-electron integrals, we need to consider one very im- 
portant result. 
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2.3 The Gaussian Product Theorem 



The Gaussian Product Theorem states that the product of two arbitrary 
angular momentum Gaussian functions on centers A and B can be written 

as 

G\G 2 = Gi(r,ai, A, /i, mi, ni)G 2 (r, a 2 , B, / 2 , ra 2 , n 2 ) 
= exp[-aia 2 (AB) 2 /7] x 



x 



'h+h 

J2 Mh, / 2 ,PA X ,PB X ) 

. i=0 

'mi+m2 

^2 / i (m 1 ,m 2 ,PA y ,PB y )^e-^ 

'n 1 +n 2 

fk(ni,n 2 ,PA z ,PB z )z k P e-^ 

I fc=o 



X 



(2.26) 



To show this, we first define the multiplicands as 



d = G 1 {r,a h A,l h m h n 1 )=x l \y r 2 1 z%e ai 
G 2 = G 2 (r,a 2) B,/ 2) m 2 ,n 2 )=x^^. 



(2.27) 
(2.28) 



Here, = r — A, etc. For primary analysis, take the angular momentum of 
G\ and G 2 to be zero, so 



Gi = e~ a ^; 



Go = e^ r *. 



(2.29) 



These are unnormalized, but normalization can be calculated as in Section 
2.1. It would be convenient if this product could be written as a third Gaus- 
sian, i.e. G\ • G 2 = G3, or 



e -a ir \ e -a 2 r% = Ke ~ir%_ 



;2.30) 
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Expand Eqn (2.30) using the definition of r^r^rp given above. 

e - ai r\-a 2 r% = eX p[- («! + « 2 )r • T + 2(«i A + « 2 B) • T 

-aiA ■ A — a 2 B • B] (2.31) 

= ifexp[- 7 (r-r-r-P + P-P)] (2.32) 

Comparing terms, 

7 = ol\ + a 2 

7 P = (a!A + a 2 B),thus P = aiA + a2B ( 2 .33) 

7 

which leads to the conclusion that 

Ke fPP = e - aiA ' A - a2B ' B (2.34) 

K = e -a 1 AA- a2 B.B+ 7 PP 

From Eqn (2.33), we expand P • P and use that to get a final expression for 
K, 

7 P P = 7" 1 [a\A ■ A + 2«ia 2 A • B + a|B • B] (2.36) 
K = exp [-ax A • A - a 2 B • B + (a\A ■ A + 2 ai a 2 A • B + a 2 2 B ■ B) / 7 ] 
= exp [(— a\A ■ A — aia 2 A • A — ai« 2 B • B — a|B • B 
+a\A ■ A + 2aia 2 A • B + a 2 2 B ■ B)7 _1 ] 

_ e -[aia2(AB 2 )/ 7 ] (2.37) 



if we define AB = (A — B). For two s-type functions, 

-aia 2 (AB 2 )/7l exp [- 7 (r - P) 2 ] (2.38) 



e -*ir\ e -a 2 r% = exp ' 



For more general Cartesian Gaussians, ones with arbitrary angular momen- 
tum, 

G X G 2 = dX^v7vlB z A z B g~ (aia2 ( AB)2/7 l e 7r1, (2.39) 



K 
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where we've used Eqn (2.38) to take care of the product of the exponentials. 
Now, x\, x\ and the like need to be considered. 

x\x\ = {x-A x ) h {x-B x ) h (2.40) 
(x-A x ) h = \{x - P x ) + (P x - A x )] h = (x P - (PA) x ) l \ (2.41) 

Using a standard binomial expansion, 



{xp - (PA)J^ = E(^( PA )^^7T— v = E(^( PA )« 



I J 

(2.42) 



Likewise, 



(x - B x f = (x P - (PB)^ = J2(x P y(PB) l r ( '? ) . (2.43) 

Using these, we can write 44 as a summation of xp to various powers. 

h+i 2 

44 = E 4A(/i, / 2 , (PA) X , (PB) X ). (2.44) 

fc=0 

The coefficient of in the product 44 ^ s gi yen by 

A(Zi, J 2 ,PA*,PBJ = E E ( - 1 ) 7 ( J ) ( 2 - 45 ) 

i=0,/i j=0,Z 2 ^ ' ^ J ' 

Perhaps more conveniently for implementing in a computational scheme, the 
constrained double sum in the above expression for fk can be redefined as a 
single sum. 

min(fc,2Zi-fc)* / J \ / 1 \ 
q=max(-k,k-2l 2 ) ^ J \ J J 

2i = k + q 
2j = k-q 

increments of 2 
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Whence we write the full Gaussian Product Theorem as Eqn (2.26). The 
derivation of Eqn (2.46) is left to the reader. 
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3 Sij — Overlap Integrals 

3.1 Overlap of primitive s— functions on different centers 

The integral we need to evaluate here is 

J 0t(r)0 2 (r)e/r = J e'^e'^dr 
Using the Gaussian Product Theorem as it appears in Eqn (2.38) 

oo oo oo 



— oo — oo — oo 

/ \ 3/2 

s 12 = e -«i« 2 (AB) 2 /7 m 

3.2 Overlap of contracted s-functions 

Take now 0i(r) to be centered on A and 02 (r) to be centered on B 

n m 

01 (r) = iYi^a,e- a ^,0 2 (r) = X>i e ~ 



m 



// V lit p 

0t(r)0 2 (r)dr = iV 1 iV 2 ^^a z 6 7 / e^e^dr 

i 3 

Examining one term in the double sum, 

/ \ 3 / 2 
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where 7^ = cti + j3j and = Q ' A ^ /3jB . So 

re 

si2 = awX!5>* 6 



n m r -I 3/2 

p—^iPj ( AB ) 2 / 7j j 7T 



(3.9) 



3.3 Overlap of primitive arbitrary angular momentum functions 

Overlap of arbitrary-/ functions: 

S12 = J Gi(ai, A, /i,mi,ni)G 2 (a2,B,/2,m2,n 2 )c/r (3.10) 

with 7 and P defined as before. Applying the fullness of the GPT [Eqn 

(2-26)], 

S12 = exp[-a 1 a 2 (AB) 2 /7]/ a; / 2/ / z . (3.12) 



where 



4 = 



/h+h 
J2 WiM^x^xVpe-^dx (3.13) 



i=0 
h+lr 



= Mh, / 2 ,PA„PB,) / x P e~^dx (3.14) 



Noting that any odd value of % produces a zero integral, and then using Eqn 
(2.22) for J x p e~ 1Xp dx, we finally write I x as 

(h+h)/2 . _ s || . . 1/2 

4= YL f2iM,PA x ,PB x y j—^f f^j . (3.15) 

The last remaining case is the overlap of two contracted Gaussians of arbi- 
trary angular momentum is left out of the consideration. It is not difficult 
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to write simple routines to compute overlap integrals over primitive Gaus- 
sian functions of arbitrary angular momentum using Eqns (3.12) and (3.15) 
and use those routines in the evaluation of overlap integrals over contracted 
functions. 
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4 Tij— Kinetic Energy Integrals 



The kinetic energy operator is — ^V 2 , or — ^(d 2 /dx 2 + d 2 /dy 2 + d 2 /dz 2 ) in 
Cartesian coordinates. So the kinetic energy integral over general, uncon- 
tracted Gaussian functions is 



12 = J ^(r)(-iv 2 )0 2 dr 
= I x + I v + L 



where we now define I x as 



4 = ~\J x'X^e-^\^)x%y^z^e-^ (4.2) 

Now we need to determine the action of the Lagrangian (or any piece thereof) 
on a particular Gaussian function. Sequentially applying the differential op- 
erator, 

%-{x%e- ai *») = l 2 x l f l e~ a ^ - 2a 2 x l * +1 e- a2X B (4.3) 



dx 

8,8,1 

dx dx 



x\e- a ^)) = l 2 (l 2 -l)x l f 2 e- a2X ^ - 2a 2 (2/ 2 + l)4 e ^ 2x| 

+4a 2 4 +2 e _a22:| (4.4) 

+a 2 (2/ 2 + \)x\e- a ^ - 2a\x h + 2 e- a ^ (4.5) 



Clearly, this is just a sum of three Gaussian functions related to the original 
by a shift of 0, 2, or -2 in the angular momentum portion, aside from some 
constants. 
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4.1 Asymmetric form of TV, 

Simply applying the results shown in Eqn (4.5) within Eqn (4.1) gives a form 
of Tij which appears as a sum of three overlap-type integrals with various 
multiplicative constants. To display the particular overlap integrals involved 
in that sum we will use a particular notation derived from the bra and ket 
notation common in physics. Let (±n| 7 denote a Gaussian where the angular 
momentum has been increased or decreased by n in the 7 coordinate. In other 
words, 

(+21, = x l+2 y m z n e- ar2 (4.6) 
Thus, given that the overlap between two Gaussians G\ and G 2 is 

J G X G 2 = (0|0), (4.7) 

the construction (0| + 2) x denotes an overlap integral between G\ and a 
Gaussian derived from G2 by incrementing the exponent of x by 2. In this 
way, we can write the asymmetric form of the kinetic energy integral using 
Eqns (4.2) and (4.5) as 

I x = a 2 (2l 2 + 1)(0|0> - 2a 2 2 (0\ + 2) x - ^"^ (Ol " 2), (4.8) 



4.2 Symmetric form of 



Time to try a different approach. Starting with the old definition of I x , 
and integrating by parts in x, 



>(r)dxdydz 



(4.9) 



2 



2 r 



dx 



dydz 



-00 




d0t(r)d0 2 (r) 
dx dx 



dxdydz 



(4.10) 
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The first term is of course zero because both 0i(r) and d(f>2(r)/dx go to zero 
as x — > ±00. So 

d(j)i d(j) 2 



1 




dx 



-dxdydz 



(4.11) 



Recalling Eqn (4.3). 



U = 



1 

2 



l 2 x h A l - 2 ai x l X +i y^z^e a ^ 



h+i 



l 2 x l f l - 2a 2 x l * +1 



■ m ^ n ^-^ r B dxdydz. 



(4.12) 



Thus we can reduce this cumbersome notation to something a little more 
convenient to code up. 



1 



■II 



l) x + 2a 1 a 2 (+l| + i; 
-ai/ 2 (+l|-l) x -a 2 /i(-l| + l): 



(4.13) 



It is somewhat more appealing, since Ty should be a symmetric matrix, i.e. 
Tij = Tji. This is an obvious truth when Eqn (4.13) is used to calculate T, 
but is not so from the asymmetric form. Good news - both expressions give 
exactly the same result. 
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5 Vij — Nuclear Attraction Integrals 
5.1 The need for a transformation 

Since the potential energy is due to Coulomb interaction of the nuclei with 
the electron in question, the operator to deal with is ^, where rc = |r — C|. 
Thus the integral we need to evaluate is 



V c = 
v ij 



[ 0*(r)-0,(r)dr 

= [ x l XyTz^e- a ^-x^z^e-^dr (5.1) 
J r c 

Since the operator does not affect the operand (0), we can combine the two 
Gaussian functions via the Gaussian product theorem, and make the final 
statement 

l m n 

•/ n (m,n 2 ,PA z ,PB z ) f x l P y^z n P e-^—dv (5.2) 

J r c 

2 

where K = e _aia2 ^ AB ^\ This is still intractable. Indeed, what do we have 
under the integral sign? 4 The usual Gaussian-like term x p ypZpe~ irp is not 
the bottleneck, it is easily representable as a product of three terms each 
depending on x, y, and z respectively. The term that prevents us from 
separating the three-dimensional integral over r into three one-dimensional 
integrals is — . Recall that 



rc 



= \/^c + yc + z c (5-3) 



At this point, we want to apply some sort of transform to the ^ to turn it into 



some sort of an exponential which can be combined with the other Gaussians 



4 Remember, it's a three-dimensional integral, r is a vector. 
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and result in resolution of the variables. There are many possibilities - the 
most commonly used are Laplace or Fourier transforms. Let's look at the 
Laplace transform solution in detail. 



5.2 Laplace transform 



Use the standard Laplace transform, 



oo 



j e- sr2 s x ' 2 ~ l ds, 





(5.4) 



where T(x) is the standard gamma function. You can just evaluate the RHS 
to confirm this. We want the instance where A = 1, thus 



rc 



7i- ' " / e sr ^s 1/2 ds. 



o 



(5.5) 



What occurs when we use this in the context of a potential energy integral 
involving only s-functions? 



V = [ e- air ^e- a ^—dv 



= K 



rc 

[ e-rt—dr 
J r c 



00 







(5.6) 



Conveniently, the Laplace transform takes r c l into a function with the ap- 
pearance of an s-type Gaussian of orbital exponent s centered at C. A second 
application of the GPT and we can switch the order of integration, evaluating 
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the integral over s second. 



oo 

2 B J1 



v - I tor* I ****** 



oo 

- K^j*r^™j**™* (5.7) 



oo 

= Kir Jdss-^ + s)- 3 ^ 3 ™ 2 /^. (5.8) 
o 

Now making the substitution t 2 = j^^,ds = ^s 1//2 (7 + s) 3 ^ 2 dt amazingly 
cancels just about everything leaving 

v= 2K 1 }^ dt 

7 J 

o 

which can be rewritten in terms of a standard error function when P ^ C, 

X 

erf (a) = J e- t2 dt (5.10) 
o 

v = ^pO erf ^ 1/2pC )- p * C < 5 - n ) 



or as 



IKir 

V = -, P = C (5.12) 

7 

The case of arbitrary angular momentum Gaussians is done likewise, only 
instead of the standard error function we get error function-like integrals 
referred to as an incomplete gamma function: 

l 

F m {T) = J t 2m e- Tt2 dt (5.13) 
o 
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Note that 



TT 



1/2 



erf(v / T),T > 



^o(T) = - 
^o(O) = 1 



2VT 



(5.14) 
(5.15) 



5.3 Fourier transform 

The Fourier transform approach is essentially very similar, the only difference 
is that — is represented as a t/iree-dimensional integral: 



A substitution of Eqn (5.16) into Eqn (5.2) and a reversal of the order of in- 
tegration allow the integral over r to be separated into three one-dimensional 
integrals. The resulting derivation is a little bit longer than in the case of 
Laplace transform, and the final expression looks a bit different, but both 
formulae give the same values, and that's what counts. 




(5.16) 
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6 Electron Repulsion Integrals 



The electron repulsion integral (ERI) is 

(ik\jl) = (ij\kl) = I #( ri )$(r 2 )— fa^fa^dv^ (6.1) 

J T\2 

ERIs are two-electron integrals as opposed to the one-electron integrals that 
we've dealt with so far. Nevertheless, the techniques that we have described 
in the previous sections can be applied to evaluate ERIs in closed form. We 
are going to sketch the strategy of a computation and present the final expres- 
sion for an ERI over four arbitrary angular momentum primitives centered 
on four different centers. 

1. We begin with the following general primitive ERI: 

0i(ri) = x l { A yTA z \A^{- a i r iA) 
M T i) = x l { B y^z^exp{-a 2 ri B ) 

Jc ». m c~ n c „ ~,2 



<Mr 2 ) = x 2 c c y^z^exp(-a 3 r 2C ) 
fa(r 2 ) = x^^4^exp(-a 4 r 2 2 ^) 

I = <f>i{ri)<f>j(ri) — fc (r 2 )<Mr 2 )dridr2 (6-2) 
J ' rw 

2. Combine fa with <pj and 4>k with fa using the GPT. 

3. The — factor in the integral demands some sort of integral transform to 
be applied to it, just like in the case of the nuclear attraction integrals: 

1 ' f " f /c- 2 e zk ' ri2 dk (6.3) 




ri2 2tt 2 

4. Insertion of Eqn (6.3) into the composite expression we've obtained in 
the first step yields a nine-dimensional integral which might look nastier 
than the six-dimensional ERI (6.2) we began with. However, we can 
switch the order of integration, and evaluate the integrals over six spatial 
coordinates (xi, X2, yi, etc.). 
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5. The resulting expression is a three dimensional integral over variables k x , 
k y , and k z , which can be evaluated quite easily. The resulting monstrous 
15-fold summation formula looks like this: 

lplq\lp + lq) l p J qUl ,u 2 t m p ,m q v p ,v q 

E G (f ) E E E G(,)F c (PQ 2 7, 7? /( 7 , + 7(Z )). (6.4) 

t' n p ,n q wi,W2 t" 

where 

C = Ip + Iq + rip + n q + nip + m q - 2ui - 2u 2 - 2v\ - 2v 2 - 
2wi - 2w 2 - t' - t" - t'" 

2 

Ki = exp(-aia 2 PQ hp) 
G(x) = (-ly^'fi^laJb^PA^PB^f^l^QC^QD,) 

x [pq y P + l g- 2u i- 2u 2- 2t ' ( ^ p ^ q y P +i q -2u 1 -2u 2 -t' 

{lp) Ul ' lp {lq) U2 ' lq W{lp + lq~ 2m ~ 2l/ 2 )!4- Ul - U2 - f/ 

X ui\u 2 \(l p - 2ui)\(l q - 2u 2 )\(lp + l q - 2m - 2u 2 - 2t')\(t')\ 
l 

F m (T) = J u 2m exp(-Tu 2 )du 
o 

< lp < la + k 0<lq<lc + ld 
0<Ul< lp/2 < U 2 < l q /2 

0<t' <(lp + l q - 2m - 2u 2 )/2 

etc. 

In the derivation we used the Fourier transformation of the — factor, but, 
likewise, the standard Laplace or Gaussian transforms may be used. 

What can we learn from Eqn (6.4)? 
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• We can immediately obtain a compact expression for an ERI over 4 
primitive s-type Gaussians: 

27T 5/2 KiK 2 2 

(ss\ss) = 7p7g ( 7p + 7g )i/2 F o( P Q 7 P 7«/(7 P + 7,)) (6-5) 

• Computing ERI of any but (ss|ss) type in closed form is not very prac- 
tical. Coding it up is relatively easy, but the efficiency is very poor, 
especially if you need to compute more than just one integral. 

• Eqn (6.4) may be rewritten as 

L 

/ = ^C m F m (T) (6.6) 

m=0 



where L = / a + ra a + n a + H h T is PQ 7 P 7 9 / (t p + 7 g ) and thus 

independent of l a , m a , n a , lb, etc. 5 "Angular" coefficients C m in turn 
do depend on the angular momentum indices l a , etc. It turns out that 
Eqn (6.6) is very general and serves as a starting point for many ERI 
evaluation methods. 

Having said this, let's look at more practical methods of computing ERIs. 



5 Let's refer to the exponents of x, y, and z as "angular momentum indices" 
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7 Practical methods of computing the electron repulsion integrals 



Up to this point we've used a notation system that we feel is the most suitable 
for manipulating closed form expressions for one-electron integrals. To move 
on further we have to enhance our notation to bring it in accordance with 
the one used in the literature. Following Obara and Saika 6 , we write the 
unnormalized Cartesian Gaussian function centered at R as 



where r is the coordinate vector of the electron, £ is the orbital exponent, 
and n is a set of non-negative integers. Sum of n x , n y , and n z will be denoted 
A and is the angular momentum or orbital quantum number. Hereafter n 
will be termed the angular momentum index. Henceforth, rii will refer to 
the z-th component of n, where % G {x, y, z}. Basic vector addition rules will 
apply to these vector-like triads of numbers, e.g. n + l x = {n x + 1, n y , n z }. 
A set of (A + l)(A + 2)/2 functions with the same A, C, and R but different n 
form a Cartesian shell, or just a shell. A set of integrals {(ab|cd)} over all 
possible combinations of functions a G ShellA, b G ShellB, etc. is termed a 
shell, or quartet, or class of integrals. For example, a (ps\sd) class consists 
of 3x1x1x6 = 18 integrals. 

The last comment before we dive in. We feel that it is beyond our ability to 
give here an exhaustive review of all available methods with all the details so 
that an unprepared person can read this document, understand thoroughly 
every concept, and be able to apply them in practice (read - write a com- 
puter code). Instead, we will try to give an overview of key ideas and most 
important algorithms and refer the reader to other, more informative and 
thorough sources. 

6 S. Obara, and A. Saika; J. Chem. Phys. 84, 3963 (1986). 



0(r;C,n,R) 



(x-R x r*(y-R y ) n y(z-R z ) 




(7.1) 
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7.1 Numerical integration methods 



The first group of methods evaluates ERIs by numerical integration similar 
to that used for integrating functions of one variable in college calculus. The 
idea of numerical one-dimensional integration is to approximate an integral 
by a finite sum: 

b f N 

/ f(x)dx = ^2f(x t )w t ; x l e [a; 6] (7.2) 

a 

Of course, we want to compute a six — dimensional integral, and our ap- 
proach will have to be a little bit different (why cannot we just use a formula 
analogous to Eqn (7.2)? Read the answer at the bottom of this page). 7 The 
trick is to represent an ERI as a one-dimensional integral over some function. 
Let's start with Eqn (6.6). After plugging the definition of the incomplete 
gamma function in and moving the summation sign inside the integral we 
obtain 

} L 

/ ^2C m i 2m exp(-Ti 2 )dt 

J ™ n 



I = 



q m=0 
1 



= J P L (t)exp(-Tt 2 )dt , (7.3) 
o 

where Ph(t) is an even polynomial of degree 2L with coefficients C m . 
Eqn (7.3) has the following form: 

b 

a 

where Ki(x) is a polynomial of degree /, w(x) is positive on the interval [a, b}. 
The theory of orthogonal polynomials offers a way of computing this type of 

7 Because the required six-fold summation would be computationally intractable 
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integrals exactly. Let a set of polynomials {S n (x)} be defined on [a; b] and 
orthogonal in the sense that 

b 

J S n (x)S m (x)w(x)dx = 5 mn , (7.5) 

a 

then integral (7.4) can be evaluated exactly by an n-point numerical quadra- 
ture formula: 

n 

I = Y,Ki(U)W^ (7.6) 

i=\ 

where n > I, t{ is the z-th positive zero of S n (x), and W{ is the corresponding 
weight factor computed as 

b 

Wi = J L t (x)w(x)dx , (7.7) 

a 

where Li(x) is the Lagrange polynomial: 

L ^ = (X~ti)---(X- tj-!){x ~ t l+ i) ■■■(x-t n ) ^ ^ 

(U — ti) • • • (U — U-i)(U — U+i) • • • (u — t n ) 

All the glorious details with the proofs can be found in any textbook on 
numerical methods. 

Fortunately, there exists a set of polynomials {Ri(t, T)}, known as Rys poly- 
nomials, that satisfy all necessary conditions: 

• R n (t,T) is an even polynomial of degree 2n in the variable t; 

• For any real T there exists an infinite set of such entities orthogonal to 
each other: 

l 

J R n (t, T)R m (t, T) exp(-Tt 2 )dt = 5 mn (7.9) 
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Hence, the ERI can be evaluated as 

n 

I = Y, p L{U)Wi , (7.10) 

i=\ 

where n is an integer greater than L/2, U is a positive zero of the n-th Rys 
polynomial 

R n (t h T) = 0, (7.11) 

and Wi is the corresponding weight which depends on T. Since Wi depends 
on T only, and every integral in a given class of integrals will share the 
same set of Wi, it is beneficial to use Eqn (7.10) to compute whole classes of 
integrals rather than a single integral. In fact, all practical strategies of ERI 
evaluation employ this so-called shell structure of integrals. 

What are potential advantages of using Eqn (7.10) instead of Eqn (6.6)? 
At first glance - none. To compute Pl one has to know coefficients C m , 
computing which is the main obstacle in using Eqn (6.6). Yet it turns out 
there exists a way to determine numerical values of Pl(U) without computing 
C m . The details are too lengthy to be presented here and can be found in a 
paper by Dupuis, Rys, and King 8 . The resulting expression is written as 

N 

I ^^2l x {ui)I y {ui)I z {ui)Wi , (7.12) 
i=\ 

where ui ~ U/(l — t 2 ) 1 ! 2 . Quantities Ii depend on angular momentum indices 
and can be computed numerically or recursively. 9 Most recently, Ishida 10 
have proposed several very efficient algorithms based on Eqn (7.12) and its 
variations, in which he computes Ii numerically and recursively. 

8 M. Dupuis, J. Rys, and H. F. King; J. Chem. Phys. 65, 111 (1976). 
9 J. Rys, M. Dupuis, and H. F. King; J. Comp. Chem. 4, 154 (1983). 

1() K. Ishida; J. Chem. Phys. 95, 5198 (1991); K. Ishida; Int. J. Quantum Chem. 59, 209 (1996). 
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7.2 Recursive methods 



"Recursive" implies computation of ERI from other ERIs of lower angular 
momentum. The first "recursive" method of computing ERIs was suggested 
by Boys. He applied the differential relation for Cartesian Gaussians 

^ = 2(0(r, C, n + l h R) - n 2 0(r, £ n - 1,, R) (7.13) 

to the expression for the (ss|ss) integral given in Eqn. (6.5) to obtain the 
corresponding (ps\ss) integral. Hence, multiple applications yield higher an- 
gular momentum ERIs. 

McMurchie and Davidson proposed recurrence relations for ERIs over Her- 
mite Gaussian functions 11 to compute the integrals over Cartesian Gaussians. 
Their algorithm is a significant improvement over Boys' method, but it re- 
quires a special transformation of the resulting integrals over Hermite Gaus- 
sians back into the Cartesian basis. 

Stagnant at the moment, the field was revitalized by Obara and Saika in 
1986. 12 Obara and Saika also used the differential relation for Cartesian 
Gaussians (7.13). Their approach treats one-, two-, and potentially n-electron 
integrals on equal footing. OS suggested the following recurrence relation for 

11 Hermite Gaussians are a special form of Gaussian functions which resemble the eigenfunctions of the harmonic 
oscillator. Their use is motivated by extreme simplicity of the differential relations for this type of Gaussians. For 
more details see ... 

12 Interesting to note, their paper was published in the April 1st issue of the Journal of Chemical Physics. No 
kidding. 
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ERIs over primitive Cartesian Gaussians: 13 

[a+ l,,b|cd] (m) = P^[ab|cd] (m) + VKP z [ab|cd] (m+1) 
+ g ([a - li.blcdjM - ^[a- l,,b|cd](- +1 )) 

+ ^ (ja,b - l,|cd]H - ^[a,b - U\cd}^ +1 ^ 



where 



oo 

2 r . / 



(ablcd)M = — dul — ^ (abided) (7.15) 
V^J \p + u z J 



o 



(ab|u|cd) = y dridr 2 0(ri; C a ,a, A)0(ri; Cfe,b,B)exp(-w 2 r 2 2 ) x 

0(r 2 ;Cc,c,C)0(r 2 ;Cd,d,D) , (7.16) 

and 

C = C + & (7-17) 

V = + (7-18) 

' = eft (7 ' 19) 

P = C ° A + aB (7.20) 

Ca + C& 

Q = • (7.21) 

Cc + u 

(ab|cd)( m ) is an auxiliary ERI playing a central role in OS manipulations. 
Note that (ab|cd)^ is a "true" ERI. 

13 The presented equation is only one of the four possible cases, the other three being the relations increasing 
angular momentum on centers B, C, and D. 
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Thus, Eqn (7.14) allows to compute any ERI from a set of (ss|ss) m integrals 

(00|00) (m) - F m (T) (7.22) 

with < m < A a + Afe + A c + A^. This result is of course equivalent to 
Eqn (6.6). However, painful evaluation of the coefficients in the sum (6.6) is 
replaced by a recursive application of a simple relation. 

Of course, the OS RR allows the shell structure of ERIs to be exploited. The 
strategy of computing primitive ERIs using OS RR is simple: 

1. compute all necessary F m (T); 

2. apply Eqn (7.14) repeatedly to compute the target (shell) of integrals. 

To compute contracted ERIs all possible combinations of primitives have to 
be evaluated separately using steps 1 and 2 and then contracted together 
to form the final value. Authors found a computer implementation of their 
algorithm to be vastly superior to other programs available at the time. 

In an attempt to improve the performance of the OS RR for contracted ERIs, 
Head-Gordon and Pople 14 suggested two complementary relations. The first 
one is termed the Vertical Recurrence Relation (VRR) and is a special case 
of the OS RR with b and d set to 0: 

[a+ li,0|c0] (m) = P^[a0|c0] (m) + TyP,[a0|c0] (m+1) 

+ £ ([a - I,,0|c0]<"»> - -^-[a- l^OlcO]^ 

+ 2K+^) [a0|c " 1< ' 01(m+1) - (7 - 23) 

The second relation is known as the Horizontal Recurrence Relation (HRR), 
also referred to as the Transfer Relation, and serves the purpose of "trans- 
ferring" the angular momentum form centers A and C to centers B and D 

14 M. Head-Gordon and J. A. Pople, J. Chem. Phys. 89, 5777 (1988). 
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respectively: 

(a,b + li|cd) = (a+li,b|cd) + AB z (ab|cd). (7.24) 

The prefactors in (7.23) still depend on the exponents of the Gaussian func- 
tions, and therefore cannot be applied to ERIs over contracted functions. 
Prefactors in HRR depend on the geometric variables only, which are com- 
mon for all primitive functions in a contraction. Hence, HRR can be applied 
to "contracted" ERIs. The computation strategy has to be adjusted accord- 
ingly: 

1. For each combination of primitive Gaussians [ab|cd] do: 

(a) Compute F m (T); 

(b) Apply VRR to build primitive [eO|fO], where A a < A e < A a + A&, 

A c < A/ < A c + A^; 

2. Contract all [eO|fO] to form contracted (eO|fO) integrals; 

3. Apply HRR to form contracted (ab|cd). 

Head-Gordon and Pople's algorithm (HGP) outperforms the standard OS 
method for contracted integrals by virtue of shifting part of the workload 
outside the contraction loops. It is considered to be nearly optimal for un- 
contracted high angular momentum functions. 

For contracted low-angular momentum classes there exist other recurrence 
relations which we will not consider here. Interested readers should refer to 
the following papers: 

P. M. W. Gill and J. A. Pople, Int. J. Quantum Chem. 40 753 (1991); 
S. Ten-no, Chem. Phys. Lett. 211, 3963 (1993); 
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8 Summary 

Authors gratefully acknowledge Jason Gonzales for help with proofreading 
this document and useful comments. E.F.V. also thanks Tool for developing 
creative atmosphere in the room. 
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9 Appendix 



Let us derive the explicit expression for the electron repulsion integral over 
primitive Gaussians using the Fourier transform for the ^ factor: 

0i(ri) = ^S^exp(-air^) 

^■(n) = Xi B y?B z iB Q M-a2r 2 1B ) 

0fc(r 2 ) = x l ^ c y^z^exp(-a 3 rl c ) 

0z(r 2 ) = x l 2 d D y™*z™ d D exp(-a 4 rl D ) 

I = 0i(ri)0j(ri) — 0fe(r 2 )^(r2)dridr 2 

7 r 12 
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